#include "normalization.h"

void initOPTS(OPTS* opt,int m,int n,int k){

	opt->m = m;
	opt->n = n;


	opt->registry[0] = '\0';
	opt->init = 2;
	opt->x0 = NULL;
	opt->c0 = 0.0;

	//Termination
	opt->maxIter = 10000;
	opt->tol = 0.0001;
	opt->tFlag = 5; // Max iterations is default

	//Normalization
	opt->nFlag = 0;
	opt->mu = NULL;
	opt->nu = NULL;

	// Method
	opt->method = 0;
	opt->lFlag = 0;

	// Regularization
	opt->rFlag = 0;
	opt->rsL2 = 0.0;

}

void setOPTx0(OPTS* opt,double* x0){
	if(strstr(opt->registry,"x0")== NULL && x0!= NULL)
		strcat(opt->registry,"x0");
	opt->x0 = x0;
}

void setOPTc0(OPTS* opt,double c0){
	if(strstr(opt->registry,"c0")== NULL )
		strcat(opt->registry,"c0");
	opt->c0 = c0;
}

void setOPTmu(OPTS* opt,double* mu){
	if(strstr(opt->registry,"mu")== NULL && mu!= NULL)
		strcat(opt->registry,"mu");
	opt->mu = mu;
}

void setOPTnu(OPTS* opt,double* nu){
	if(strstr(opt->registry,"nu")== NULL && nu!= NULL)
		strcat(opt->registry,"nu");
	opt->nu = nu;
}


int hasGot(OPTS *opt,char * field){
	if(strstr(opt->registry,field)== NULL)
		return 0;
	return 1;
}

void initNormalization(OPTS *opt,double *A){

	//%%%%%initialization of one vector can be moved elsewhere
	int maxval = maxof2(opt->m,opt->n),i,j;
	double* oneVector;
	double norm2Val;
	oneVector = (double*) malloc( sizeof(double)* maxval );
	for(i=0;i<maxval;i++)
		oneVector[i] = 1.0;
	if(opt->nFlag !=0){
		if(!hasGot(opt,"mu")){
			opt->mu = (double *) malloc(sizeof(double)* opt->n);
			dgemv('T',opt->m,opt->n, (1.0/((double)opt->m)), A,opt->m,oneVector,1.0,0.0,opt->mu,1);
		}
		if(opt->nFlag ==1){
			if(!hasGot(opt,"nu")){
				opt->nu = (double *) malloc(sizeof(double)* opt->n);
				for(i=0;i<opt->n;i++){
					norm2Val = dnrm2(opt->m,(A+i*opt->m),1);
					(opt->nu)[i] = sqrt( norm2Val * norm2Val  / ((double)opt->m));
				}
				if( (opt->nu)[i] < 1.0e-10 )
					(opt->nu)[i] = 1.0;
			}
		}
		else
			if(!hasGot(opt,"nu")){
				opt->nu = (double *) malloc(sizeof(double)* opt->m);
				for(i=0;i<opt->m;i++){
					(opt->nu)[i]=0.0;
					for(j=0;j<opt->n;j++)
						(opt->nu)[i] += A[i +j* opt->m]*A[i +j* opt->m];
					(opt->nu)[i] = sqrt((opt->nu)[i]/ ((double)opt->n) );
					if( (opt->nu)[i] < 1.0e-10 )
						(opt->nu)[i] = 1.0;
				}
			}
	}
}


void normalizedmv(char trans,OPTS* opt,double* A, double* x, double* y){
	int maxval = maxof2(opt->m,opt->n),i,j;
	double* oneVector;
	double sumx;
	double * invNu;
	oneVector = (double*) malloc( sizeof(double)* maxval );
	if(opt->nFlag == 0)
		dgemv(trans,opt->m,opt->n,1.0,A,opt->m,x,1,0.0,y,1);
//elseif (opts.nFlag==1)
//    ATy=A'*y - sum(y) * mu';  ATy=ATy./nu;
//else
//    invNu=y./nu;              ATy=A'*invNu-sum(invNu)*mu';
//end
	
		if(trans == 'T' || trans == 't'){
			if(opt->nFlag ==1){
				sumx = ddot(opt->m,x,1,oneVector,1);
				dcopy(opt->n,opt->mu,1,y,1);
				dgemv(trans,opt->m,opt->n,1.0,A,opt->m,x,1,(0.0-sumx),y,1);
				for(i = 0;i<opt->n; i++)
					y[i] = y[i] / (opt->nu)[i];
			}
			else{
				invNu = (double*) malloc( sizeof(double)* opt->m );
				for(i = 0;i<opt->m; i++)
					invNu[i] = x[i] / (opt->nu)[i];
				sumx =  ddot(opt->m,invNu,1,oneVector,1);
				dcopy(opt->n,opt->mu,1,y,1);
				dgemv(trans,opt->m,opt->n,1.0,A,opt->m,invNu,1,(0.0-sumx),y,1);
			}
		}
//elseif (opts.nFlag==1)
//    invNu=x./nu; mu_invNu=mu * invNu;
//    Ax=A*invNu -repmat(mu_invNu, m, 1);
//else
//    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
//end
		else{
			if(opt->nFlag ==1){
				invNu = (double*) malloc( sizeof(double)* opt->n );
				for(i = 0;i<opt->n; i++)
					invNu[i] = x[i] / (opt->nu)[i];
				sumx =  ddot(opt->n,invNu,1,opt->mu,1);
				dcopy(opt->m,oneVector,1,y,1);
				dgemv(trans,opt->m,opt->n,1.0,A,opt->m,invNu,1,(0.0-sumx),y,1);
			}
			else{
				sumx = ddot(opt->n,x,1,opt->mu,1);
				dcopy(opt->n,oneVector,1,y,1);
				dgemv(trans,opt->m,opt->n,1.0,A,opt->m,x,1,(0.0-sumx),y,1);
			}
		}
}